Phase transitions on Markovian bipartite graphs — an application of the zero-range 

process 



Otto Pulkkinen* and Juha Merikoskit 
Department of Physics, P.O. Box 35 (YFL), FIN-40014 University of Jyvaskyla, Finland 

(Dated: February 2, 2008) 

We analyze the existence and the size of the giant component in the stationary state of a Markovian 
model for bipartite multigraphs, in which the movement of the edge ends on one set of vertices of 
the bipartite graph is a zero-range process, the degrees being static on the other set. The analysis is 
based on approximations by independent variables and on the results of Molloy and Reed for graphs 
with prescribed degree sequences. The possible types of phase diagrams are identified by studying 
the behavior below the zero-range condensation point. As a specific example, we consider the so- 
called Evans interaction. In particular, we examine the values of a critical exponent, describing the 
growth of the giant component as the value of the dilution parameter controlling the connectivity 
is increased above the critical threshold. Rigorous analysis spans a large portion of the parameter 
space of the model exactly at the point of zero-range condensation. These results, supplemented 
with conjectures supported by Monte Carlo simulations, suggest that the phenomenological Landau 
theory for percolation on graphs is not broken by the fluctuations. 
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I. INTRODUCTION 

The common way of implementing time dependence 
in the models of random graphs is to couple the time 
and the number of edges in such a way that, starting 
with an empty graph with a given number of vertices, a 
new edge is added uniformly at random to one of the 
vacant edge locations at each time step. This proce- 
dure is known simply as the graph process [1]. During 
the past few years, a large variety of models generalizing 
the graph process have appeared in the physics literature 
(see Ref. [2] for a review), most of them motivated by 
the structure of the Internet and the World Wide Web. 
One of the central ideas of these models is known, by 
the work of Barabasi and Albert [3], as preferential at- 
tachment. In the Barabasi- Albert model new vertices 
are introduced one at a time and joined to one of the 
existing vertices chosen with probability proportional to 
the degree of the vertex, i.e. the number of edges already 
joined to it. This preferential attachment leads to what 
Price [4] calls cumulative advantage, and it can be used to 
construct so-called scale-free networks, i.e. graphs with 
power-law distributed degrees. 

In the graph process and the Barabasi- Albert model, 
the edges, once attached to the vertices, retain their po- 
sitions ad infinitum. In the present article, we consider 
what happens if internal restructuring by rewiring of the 
graph with the existing edges dominates the dynamics, 
that is, the time scales associated with addition and re- 
moval of edges and vertices, have become very long as 
compared to the times between the rewiring events. The 
dynamics is taken to be the simplest possible capable 
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of providing the cumulative advantage - namely of zero- 
range [5] so that the rate of a rewiring event, i.e. a jump 
of an edge end from one vertex to another, depends only 
on the degree of either the initial or the final vertex of 
the move. We focus on the connectivity properties of 
the graphs in the stationary state. The question about 
the existence and size of the giant component, a con- 
nected cluster occupying a macroscopic fraction of ver- 
tices, turns out particularly interesting because of the 
condensation phenomena induced by the zero-range pro- 
cesses. Closely related models have been studied on a 
phenomenological level by Palla et al. [6] and rigorous 
results for models of polymerization sharing the same 
kind of ideas have been obtained by Pittel et al. [7, 8]. A 
statistical mechanics approach to a rather general class of 
reversible graph- valued processes is presented in Ref. [9] . 

We take bipartite multigraphs as the state space of the 
problem so that the system consists of external agents 
with static degrees and another set of vertices, for exam- 
ple collaboration events, on which the movement of the 
edge ends, coming from the set of agents, takes place. 
A process with unipartite states is obtained as a special 
case to our model, but its behavior is not as rich as in 
the bipartite version. 

The structure of the article is as follows: In the fol- 
lowing section, we give a precise definition of the model 
and discuss the reasons for choosing the particular state 
space. Section III reviews some fundamental results 
about zero-range processes and a Poisson approxima- 
tion for the numbers of vertices of given degree is con- 
structed. This is applied to subcritical zero-range pro- 
cesses in section IV. In section IV A, two simple exam- 
ples are given and general implications of the results for 
subcritical zero-range processes are discussed in section 
IV B. In section V, we concentrate on restructuring pro- 
cesses with the so-called Evans interaction. A detailed 
view of the phase diagram and the critical exponents is 
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FIG. 1: (a) Example of a bipartite multigraph with vertex 
sets V and W. The static degrees of v £ V are denoted by 
K v and the evolving degrees of w £ W by X w — X w + X w 
with X w and X w the numbers of edges leading to vertices of 
degree one and two, respectively. The arrow and the dotted 
line show a possible transition, (b) The projection of the 
bipartite graph in (a) to the vertex set W . 

given in subsection VA. The connection with a model 
of diluted scale-free networks is also established. Con- 
clusions are made and some open problems presented in 
section VI. 



II. STATEMENT OF THE PROBLEM 

We shall consider bipartite multigraphs with vertex 
sets W and V such that \W\ = L and \V\ = M 1 + M 2 
(see Fig. 1(a)). Here Mi denotes the number of vertices 
of degree 1 and M 2 the number of vertices of degree 2 and 
these degrees are assigned to randomly chosen vertices in 
V independently of each other. This set up on V will 
be static during the time evolution. The total number of 
edges is N = Mi + 2M 2 and we shall consider only cases 
with the densities 

r :=2M 2 /Ne [0,1], p := N/ L e (0, oo) (1) 

fixed on the passage to the limit L — > oo. 

Let us then discuss the dynamics. Initially the N edge 
ends are chosen to be distributed uniformly and indepen- 
dently of each other on the vertices in W. The restructur- 
ing of the graph takes place in continuous time by jumps 
of the edge ends on this set. We take the process on W to 



be a zero-range process [5, 10] with Mi particles of type 
1 (edges to vertices of degree 1 in V) and 2M 2 particles 
of type 2 (edges to vertices of degree 2 in V), so that a 
vertex w £W loses a randomly chosen edge end after an 
exponentially distributed time, where the parameter of 
the exponential distribution is a function of the degree 
of w. More precisely, if w € W is an end vertex of X^ 
edges to vertices of degree 1 and of X%, edges to vertices 
of degree 2, the exponential rate is given by g(X^ + X^), 
where g is a given positive function with bounded incre- 
ments, and the edge end to jump is one of those leading 
to vertices of degree 1 with probability X^/(X^ + X^) 
and one of those leading to vertices of degree 2 other- 
wise. Immediately after leaving the vertex w, the edge 
end joins one of the other vertices according to a symmet- 
ric, irreducible transition matrix (P WlZ ), w,z £ W. The 
stationary distribution of this zero-range process can be 
found explicitly and will be studied in the next section. 
Clearly, the stationary distribution of the corresponding 
bipartite multigraph -valued process is such that all the 
states with the same degree sequence are equiprobable. 

In this article, we study the structure of the graphs 
in the stationary state of the process. In particular, wc 
concentrate on the properties of the projection of the 
bipartite graph, also called one-mode network [11], on the 
set W, obtained by adding an edge between two vertices 
in W if one can be reached from the other by traversing 
two edges on the original bipartite graph, see Fig. 1(b). 
In other words, the projection tells us how the vertices 
in W are mutually connected by the vertices of the other 
type. The particular questions we ask are then, when 
does a giant component, a connected subgraph of size 
proportional to L, exist on the projection and precisely 
of how many vertices such a component is made? 

There are several reasons for making the particular 
choices concerning the state space. First of all, the 
zero-range dynamics is easily implemented in this set- 
ting and it contains the unipartitc model as the special 
case Mi = 0. The bipartite structure really brings some- 
thing more to the theory: The Mi vertices of degree 1 
are dangling ends from the point of view of connectivity, 
and do not appear on the projection graph, but they do 
have influence on the jump rates of the zero-range pro- 
cess. The parameter r = 2M 2 /N controls the amount of 
this 'dark matter' and lowering its value while keeping 
p = N/L fixed corresponds to diluting the set of edges 
on the projection. Therefore it will be referred to as a 
dilution parameter. Varying the values of r and p si- 
multaneously in such a way that their product remains 
constant can be used to sweep the phase diagram of the 
underlying particle system, characterized by the density 
p alone, without changing the density of edges on the 
projection. In addition, the restriction of the static de- 
grees to at most 2 allows us to use the usual configuration 
model [1], and therefore the results of Molloy and Reed 
[12, 13], directly without going through the proofs in the 
bipartite setting. Moreover, we believe that the extension 
that allows higher but bounded degrees or even degree 
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distributions with exponential tails does not change the 
qualitative features of the phase diagrams or the values 
of the critical exponents. 



III. ZERO-RANGE PROCESSES AND POISSON 
APPROXIMATIONS 

Let the dynamics of the edge ends on the set W be 
as described in the previous section and let ^Il,m u m 2 
be the set of ordered partitions of two numbers Mi and 
2M 2 both in L non-negative parts. Then, for (77, £) G 
^l,Mi,m 2 , one shows by verifying that the detailed bal- 
ance condition holds that the stationary (canonical) dis- 
tribution of the zero-range process is given by 



Mi,Mi,M 2 (??,0 

1 



(2) 



n 



Z^(L,M 1 ,M 2 ) 11 V 



where g\(n) = YYk=i 9$) ^ s a generalized factorial with 
the convention that g\(0) = 1 and Z^(L,Mi,M 2 ) is a 
normalization factor (canonical partition function). Let 
D m ,n denote the number of vertices in W with exactly 
to edges leading to vertices of degree 1 and n edges lead- 
ing to vertices of degree 2. Then (2) implies an image 
measure 



C r L,Mi,M 2 (^) = 



1 



M 1 M 2 

><nn 

m=0 n— 



Z a (L,M 1 ,M 2 



m + n 



(3) 



d n 



m ) g\(m + n) 



for the joint random variable D — (D m ^ n ) in the station- 
ary regime. In Eq. (3), d is such that the three condi- 

tl0nS T,m=0 T,n=0 d ™,n = L, Y, m =Q T,n=o( m + n ) d ™,n = 

Mi + 2M 2 and J2mLo Y.n=o m dm,n = Mi are satisfied. 
Otherwise <tl,m 1 ,m 2 vanishes. 

Our purpose is to get a grip of the random variables 
D n := J2m=o D m ,n, i-e. the numbers of vertices of degree 
n on the projection to W. More formally V — (D n ) 
is the asymptotic degree sequence of Molloy and Reed 
[12, 13] in our problem and we have to show that, with 
high probability [14], we get an asymptotic degree se- 
quence T>, which is in their sense well-behaved (see ap- 
pendix A for definitions and discussion). Because of the 
constraints on d in Eq. (3), the random variables Drn,n 
are not independent, and it would be difficult to extract 
the behavior of D n from the measure ol,M\Mi directly. 
Next we are going to construct a sequence of independent 
variables that is shown to be a good approximation for 
the original sequence. In other words, we shall study the 
problem in the grand canonical ensemble and prove that 
the canonical and grand canonical ensembles are equiva- 
lent. 

For zero-range processes, approximations by indepen- 
dent variables have been considered by Grofikinsky et 



al. [15] and Jeon et al. [16]. In their highly influential ar- 
ticle, Jeon et al. were able to give a rigorous proof for the 
behavior of the largest cluster not only for subcritical pro- 
cesses and for the cases with b > 3 in what is sometimes 
called Evans interaction [10], g(k) = l+b/k+0{kr { - 1+s ^), 
but also for a large class of vanishing rates. Within this 
class, which essentially consists of functions that vanish 
faster than exp(— tVlog k), a single vertex would hold 
all but of the order Lg(L) edge ends in our model, so 
that the projection of the graph would be in a flower- 
like state with one vertex having many self-loops. We 
choose to study the vanishing cases no further. By re- 
stricting ourselves to the cases with the interaction func- 
tion g bounded away from zero, we also avoid the trunca- 
tion of the series involved in the grand canonical analysis. 
The phase diagram of the Evans model was also discussed 
by Grofikinsky et al. and heuristic arguments, supported 
by simulation results, for the dynamics of condensation 
were given. The results of these articles concerning the 
Evans interaction will be applied in later sections. The 
general method and the validity of approximations by in- 
dependent variables is discussed in a delightful paper of 
Arratia and Tavare [17]. 

Let now g be bounded away from zero and let 
C m .n, to, n € Z+, be independent Poisson(A m ,„(a, 7, <j>))- 
distributed random variables, where 



A m ,n(a,7, 4>) 



1 /m + n\ a(l -7) m 7 n </> m +™ 
~Z{4>)\ to J g\{m + n) ' ( ' 



m = E 



k=0 



gm : 



(5) 



and a, 7 and <fi are real parameters whose values will be 
determined later. Here Z is the grand canonical partition 
function. Its radius of convergence, which we shall denote 
by is positive because g is bounded away from zero. 
In order to avoid cumbersome expressions, we set 



00 00 



Ai 

A 2 
A 3 



m—0 n—0 

00 00 

m— n—0 

00 00 

{ zZzZ mCm ' n = Ml }- 



(6) 



m— n—0 



Then one can show by a direct calculation that 
(Ad,o, • • -,D Ml ,2M 2 ) = (C ,o, • • • \AinA 2 nA 3 ) in law and, 
because of cancellation in the conditional probabilities, 
independent of the values of a, 7 and 4>. So (C TOi „) really 
is a candidate for a good approximation. The construc- 
tion suggests the following for the approximative system: 

1. A direct calculation shows that the system size 
Dm=o Y^=o Cm.n is a Poisson(a)-distributed ran- 
dom variable. 
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2. Again by a direct computation, for t £ 



Ecxp it J2( 



m— n— 



in that 



^ y^(m + n)C„ 



Vm-0 n— 



in law, with the variables Y w independent and dis- 
tributed according to the grand canonical distribu- 
tion 



v<t>{k) :-- 



1 



Z(<t>) g\{kY 



k £ Z+. 



(7) 



In other words, v$ controls the vertex degrees lo- 
cally. The expectation of a ^-distributed variable 
will be denoted by R(4>): 



R{4>) 



EY W 



4>Z'(<I>)/Z(<I>). 



(8) 



The parameter (f> is known as the fugacity. 
3. In a similar manner one shows that 



E 



exp ( it mC n 



A, 



so that, given A\ and A 2 , the number of edges lead- 
ing to vertices of degree equal to unity has the bi- 
nomial distribution with parameters N and 7. 

Now suppose that one would like to get an exponen- 
tially decaying upper bound for the probability of an 
event B concerning the dependent variables {D m , n ), and 
let B be the same event for the independent variables 
(C m , n )- Suppose further that the parameters a, 7 and 
4> can be chosen in such a way that the probability of 
the conditioning event A\ n A 2 fl A3 is not exponen- 
tially small. Then it suffices to get the upper bound 
for the independent case: P(B) = P(B\A 1 C\A 2 C\ A 3 ) < 
P(B)/P(A 1 n A 2 n A 3 ). The use of this conditioning 
device is nicely illustrated in the articles of Corteel et 
al. [18] and Jeon et al. [16]. 

How to best choose the values of the parameters? The 
point 1 in the list above tells that the approximative sys- 
tem size X)m=o S^Lo Cm," i s sharply peaked at a with 
a standard deviation of \fa so that the simplest choice 
is to match the parameter with the size of the canonical 
system L. Notice that the choice is not unique because 
terms of the order o(L) could be added without affecting 
the behavior in the large system size limit. The points 2 
and 3 in the list show that, with the appropriate condi- 
tionings, the same type of argument applies to the two 



other random variables in the conditions (6) as well. We 
thus set 



a = L, 7 = 2M 2 /N = r 



and, if possible, 



R((j)) = N/L = p. 



(9) 



(10) 



Clearly, the first two equations can always be fulfilled but 
in the third one the case 



lirn R(<j>) < p 



(11) 



marks a condensation transition. Roughly speaking, the 
best one can do in that case is to set (j) = 4>, so that 
locally one observes degrees distributed according to z/$, 
just like in a system with density p c . Then the rest of the 
mass, that is (p — p c )L edges, must hide in the o(L)-sets 
of vertices since those edge ends are not picked up by the 
approximating measure. It was proven in Refs. [15, 16] 
that this really is the true picture and, furthermore, the 
authors of Ref. [15] claim that the set of o(L) vertices 
with the rest of the mass is in fact a single vertex for a 
large class of rate functions. Especially, the Evans model, 
which we shall study in section V, belongs to this class. 
As a by-product of the discussion above, we get some- 
thing really important to the theory. From the point 
of view of our application the condensates for vanishing 
rates and the ones bounded away from zero are com- 
pletely different: In the first case, the components on 
the flower-like projection graph are very small, while in 
the second, the condensate on a single vertex implies the 
existence of a giant component. 

Recall that D n := J2m=o D m ,n is the number of ver- 
tices of degree n on the projection, so that C' n :— 
12m=o Cm.n is an independent approximation for this 
quantity. We define 

which equals the expected value of C n divided by the 
system size. To show that the conditions of Molloy and 
Reed are satisfied, we first prove that the fluctuations C n 
around v^(n)L arc not too large. The following inequal- 
ity is rather easily obtained for any e > using Chcrnoff 's 
bounds (see e.g. [19]): 



> e 



< exp (-eL/(e/!/J(n))) , (13) 



where f(x) = ((1 + x) log(l + x) — x) /x, with x > 0, 
is strictly increasing. Using here the trivial bound 
/(e/fjW) > /(e) and the fact that e/(e) > (l + e)e/(l + 
e/2) — e = e 2 /(2 + e) yields that, with high probability 
for any S > 0, 



C n = Mn)L + o(L 1 / 2+s ). 



(14) 



In order to get the estimates for the original sequence, 
one has to bound P{A\ fl^fl As) from below. 
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IV. SUBCRITICAL ZERO-RANGE PROCESSES 

Let us now derive a subexponential lower bound for the 
probability of the conditioning event in the subcritical 
cases p c := lim^^j, <j>Z'((f>)/Z((j>) > p. Notice first that 
with the choices made in Eqs. (9) and using the results 
for the variables in Ai and A3, the probability of A\ 
is bounded below by c\j\[L for some constant c\ and 
also P(A 3 |Ai n A 2 ) > czj\fL. So we have to bound the 
probability of A 2 given A\. By subcriticality and from 
the fact that 4>Z' '((/)) / 'Z((f>) is strictly increasing it follows 
that there is a unique 4> < $ for every value of the density 
such that Eq. (10) is satisfied and the distribution u$ 
therefore has an exponential tail. Therefore the usual 
local limit theorem [20] applies and the probability is 
at least c 2 /VX. Thus P(A X n A 2 n A3) > cL~ 3 / 2 and 
the bound (14) for the fluctuations of C n applies to D n 
as well. Moreover, the tail of i/J is exponentially thin 
because v r ^(n) = 0(V(j,(n)/r) and v§ shows exponential 
decay. Also, since J2k>n l '4>(^) = 0{c n ), for some c € 
(0,1) and large n, the largest degree is O(logL) with 
high probability. Now it is easy to check that, with high 
probability, the degree sequence (D n ) is such that the 
conditions of Molloy and Reed, reproduced in appendix 
A, are satisfied. 

Next we give the results for the subcritical cases. Let 
$ > and p c > p and define 



<j>z"{<j>y 



(15) 



Then, the parts of Theorem 1 of Ref. [12] concerning 
the existence of the giant component and Theorem 1 of 
Ref. [13] state that, with high probability: 

1. If r < r c ((/>), the largest component on the projec- 
tion to W has at most of the order log 3 L vertices. 

2. If r > r c ((f>), there is exactly one component on 
the projection with more than TlogL vertices for 
some constant T, and the size of that component 
is LA + o(L), where 



A = 1 



Z((l - W) 

m 

and (3 is the positive solution to 
Z'((l - /?r)0) 



0=1- 



(16) 



(17) 



Notice that, since Z'({1 — /3r)<f>) is a convex function 
of (3 on (0,1), Eq. (17) has a unique positive solution 
exactly when r > r c ((f>). 

We can also deduce the leading order of A as a function 
of r — r c ((f>): The expansion of the right-hand side of 
Eq. (17) for small (3 (all the derivatives of Z(<f>) exist 
because < 3>) yields 



0> 



r 2 (t) 2 Z"'((j)) \r c (<j>) 



1 



(18) 



as r — > r c from above, so that now expanding the right- 
hand side of Eq. (16) to the first order in (3 we have 



A w rp(3 • 



i<pz"{<py 

Z(<t>)Z>»(<j>) 



(r - r c {4>)). 



(19) 



Thus the growth starts linearly. Furthermore, the 
uniqueness of the solution to Eqs. (16) and (17) allows 
us to calculate the leading order of A when the critical 
curve is crossed by increasing the density p with constant 
r. At the point (<fi + e, r c {4>)) with <fi + e < $, we get 



wz"{4>y 

z{4>)Z'"{4>) 



(20) 



for the size of the giant component. Note that, in case of 
a finite radius of convergence <& of the partition function, 
the 4> $ limits of the amplitudes in expressions (19) 
and (20) depend crucially on the existence of the third 
moment of a ^-distributed random variable. 



A. Two examples 

As a first example, let us study the case of non- 
interacting random walks. Then g(k) = k and Z(<p) = 
exp(0), so that from Eqs. (8,10) and from the definition 
of the critical density of the zero-range process, Eq. (11), 
we get that <j> = p and p c = 00. Also, from Eq. (15), it 
follows that r c (p) = 1/p, for p > 1, and there is no phase 
transition for p < 1 (remember that r = 2M 2 /N 6 [0, 1]). 
Of course, by r c (p) we mean r c ((f>(p)) with <f>(p) the solu- 
tion to Eq. (10). Setting c := rp one recovers the results 
familiar from the model G n ^M=cn of random graphs [1]: 
The phase transition occurs at c = 1 and the size of the 
giant component is LA + o(L) with A given by the solu- 
tion to A — 1 — exp(— cA). Furthermore, Eq. (19) is con- 
sistent with the known expansion A = 2e — 8e 2 /3 + 0(e 3 ) 
for c = 1 + e. 

In our second example we take the jump rates to be 
degree independent, g(k) = 1. Now the partition func- 
tion has a finite radius of convergence, $ = 1, but still 
there is no condensation transition, so that p c = 00. In 
this case we have r c (p) = l/(2p), for p > 1/2, and no 
transition for smaller densities. Quite surprisingly, the 
size of the giant component has a simple expression 



4 

rp ' 



from which, again in harmony with Eq. (19), we get A = 
4pe/3 - 56(pe) 2 /27 + (9(e 3 ), for r = l/(2p) + e. 



B. General remarks on the phase diagrams 

The critical curve can be alternatively written in terms 
of the function R(<f>) defined in Eq. (8): 



r c (0) = R{4>) - 1 + 



R 



IV 1 



R(4>)J 



(21) 
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from which it follows that 



P-1 



(22) 



The fact that the equation for the critical curve involves 
the second derivative of the partition function, or equiva- 
lcntly the first derivative of R, already tells that the phase 
diagram of the graph-valued process can be more com- 
plicated than for the zero-range process itself, in which 
the existence of the phase transition is simply dictated by 
the asymptotics of the ratio Z' \4>) / Z (cf>) . Three different 
kinds of behavior can be extracted from the results (15) 
and (21): 

1. When R((f>) — > oo as <j> — ► <&, i.c p c is infinite, the 
critical curve decreases to zero as a function p in 
such a way that r c (p) > for p arbitrarily large. 
The high density asymptotic behavior of the curve 
is either r c (p) <~ 1/p or given by the nontrivial part 
in the denominator of (21). This point will be ex- 
plored further shortly. 

2. When the first derivative of the partition function 
converges as <f> — > $, in that p c is finite, but the 
second derivative diverges, we have a condensation 
transition, and the critical curve approaches zero 
as p tends to p c from below. Therefore, a giant 
component exists at and above the zero-range crit- 
ical point for any non-zero r. Near the condensa- 
tion transition, where <f> is just below <&, the critical 
curve goes like 



(23) 



3. If also the second derivative of the partition func- 
tion converges, lim p _>p c r c (p) is non-zero. In fact, 
the values of the dilution parameter r being re- 
stricted to the unit interval, this case is further 
divided into two categories according to whether 
this non-zero limiting value is greater or less than 
1. For lim p ^ Pc r c (p) > 1, no phase transition can 
be achieved by varying the dilution parameter r. 
In both cases, however, if the condensation occurs 
on a single vertex, we know that above p c there 
is a giant component, but the theory developed so 
far tells nothing about its existence exactly at the 
zero-range critical point. 

These three points illustrate the shape of the critical 
curve at the high density limit. We remark that, for 
small values of p, the curve is expected to pick up the fea- 
tures of the non-interacting case g(k) — k, provided that 
lim p ^ Pc r c (p) < 1. An example of a concrete phase dia- 
gram, with all the features discussed here, will be given 
in section V. 

We saw that the interactions in the two examples of 
section IV A produced critical curves inversely propor- 
tional to the density. Let us now deduce, for which func- 
tions g this proportionality holds, at least for densities 



high enough. For this purpose, set 

rc{<t>) = J(fy 



(24) 



with a > 0, in (21). This yields a linearizable first or- 
der differential equation for R, with the initial condition 
R'(0) = l/g(l), the solution to which is 



R{<j>) 



4> 



S(l) + (l-a)0' 



(25) 



where (f> € [0,g(l)/(a — 1)) for a > 1 and € [0, oo) 
otherwise. Furthermore, solving for Z(<f>) from Eq. (8), 
we get 



z{4>) 



i + 



1-a 

W) 1 



(26) 



which has the expansion 

av ' m>2 vav ;/ k=2 

(27) 



k=2 



so that 



g{k) 



.9(1), 

g(i)fc 

_ (fc-l)o-(fe-2) ' 



for k = 1 
for k > 2. 



(28) 



This form of interaction function covers a range of mod- 
els, which we now classify for different values of a: 

For a = 1 or 2 we get the examples of section IV A. 

For 1 < a 7^ 2, the tail of g(k), which determines the 
high density behavior, can be expanded in powers of 1/fc: 



\™ (a - 1)« 



(29) 



71=1 
OO 

E 

n— n*+l 



where n* = max{0, [^5f]}, and we have chosen g(l) = 
a — 1 to make the radius of convergence equal to unity. 
To leading order in the large k limit, g(k) is then given 

by 



g(k) 



l-b(a)/k, for a e (1,2), 
l + b(a)/k, for a > 2, 



(30) 



where 6(a) = |(a - 2)/(a - 1)|. 

In the first case of 1 < a < 2, which interpolates be- 
tween the cases of non-interacting edges and the degree 
independent rates, the range of 6 as a function of a is R+ , 
so that the behavior 



r c {p) 



1 



(2 + b)p 



(31) 



is expected to hold for all b <E (0, oo) in the high density 
limit for the interaction g(k) = 1 - b/k + 0(fc"( 1+(5 )). 
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The interesting case of a > 2 with decreasing interac- 
tion function is in [15] referred to as the Evans interac- 
tion due to his work presented in [10]. Now the image 
of {a > 2} under b is the interval (0, 1) — just a tiny part 
of the nontrivial parameter space < b < oo. Thus, we 
predict 



r c (p) = 



l-b 

(2 - b)p 



(32) 



for high densities in the Evans model with < b < 1 and 
something else for the rest of the phase diagram. This 
model is known to exhibit a condensation transition for 
b > 2, in which case our critical curve r c (p) will hit zero at 
p c . The region 1 < b < 2 is therefore special - the critical 
curve is strictly above zero, but falls of faster than the 
first inverse power in p. This seems to be consistent with 
the observation of Grofikinsky et al. [15] that a kind of 
precursor of condensation transition is present on this 
interval of the parameter space. The Evans model will 
be the topic of the next section. 

In the region < a < 1, the cases with a € 
((0, 1) \ {1 — l/(n+ 1) :n£ N}) can be ruled out imme- 
diately, since they lead to interaction functions breaking 
the positivity assumption. So we are here left with the 
set a e {1 — l/(n + 1) : n € N}. But this time the as- 
sumption on bounded increments is violated: g(k) would 
be finite for k € {1, . . . ,K}, where K = 1/(1 — a), and 
infinite for larger values of k. Thus we don't expect to 
find l/(ap)-behavior with a < 1 in our model. However, 
it is possible to construct processes that have station- 
ary distributions corresponding to the partition functions 
(26) by altering the dynamics only slightly. This leads to 
the notion of generalized exclusion processes [21] and, in 
these processes, the jumps to sites already occupied by K 
particles are simply suppressed. Especially, the version 
with K = 1 is the simple exclusion process [22-24]. We 
would like to remark that, since for generalized exclu- 
sions the condition p < K must be satisfied, the bound 
(22) would not be broken. 



V. EVANS INTERACTION 

In this section we take g(k) — 1+b/k, b > 0, so that the 
partition function is given by a hypergeometric function 
[25], 



m = ±(^±^.^ =F(W + bM , (33 ) 



T(k + l + b) k\ 



fe=0 



with the radius of convergence $ = 1. Also, from the 
definition (8), 



R{<P) 



4> F(2,2;2 + b;< 
~l + b ' F(l, !;! + &;< 



(34) 



from which the value of the critical density can be calcu- 
lated [15]: 



oo, for b < 2 

1/(6-2), for b > 2. 



(35) 



The formula for the critical curve is, by the definition 
(15), 



2 + b F{2,2-2 + b- 



F(3,3;3 + 6; 



(36) 



Let the size of the largest cluster in a zero-range pro- 
cess on L sites be denoted by Z* L . The following facts for 
Evans interaction are taken from Ref. [16]: 

1. When p < p c , that is tfi < 1, the distribution v$ has 
an exponential tail, and Z* L is at most of the order 
log L with high probability. 

2. At the critical point, that is <j> = 1 and b > 2, 
the distribution has a heavy tail, vi(k) ~ k~ b , 
and Z* L is exactly of the order L 1 /( b_1 ) with high 
probability. Furthermore, 



P(A 2 | Ai) > cL 1 - 



(37) 



3. For p > p c and b > 3 (this is needed in the proof), 



(38) 



in probability as L —* oo , and the size of the second 
largest cluster is o(v / L log 2 L) with high probabil- 
ity. 

The results of the previous sections cover the subcrit- 
ical cases p < p c only. As already discussed, the case 
number 3 in the previous paragraph implies the existence 
of a giant component, but its size is not known. The only 
cases that the existence of a macroscopic component is 
not yet proven are the ones exactly at the zero-range 
criticality with b > 3. This is because the second mo- 
ment calculated from the distribution v\ exists in this 
region only, and the critical curve r c (<f>) was inversely 
proportional to the second derivative of Z(4>). Therefore 
lim^! r c {4>) = if and only if b < 3. 

Next we show that, at the zero-range critical point 
with b > 3.51, our degree sequence is with high proba- 
bility such that it satisfies the conditions of Molloy and 
Reed. Evidently, the construction of the lower bound 
for the probability of the conditioning event follows the 
lines of the subcritical analysis. The result of Jeon et 
al. in Eq. (37) gives an estimate for P(A 2 \ Ai) in the 
whole phase diagram with condensates, that is b > 2, 
but due to the restriction to the region with existing sec- 
ond moment of the measure v\, things are settled by 
the usual local limit theorem. Concerning the upper 
bounds for the deviations of n(n — 2)C n /L around its 
mean n(n — 2)v{{n) (see appendix A), we see that locally 
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the condition b > 3.01 would suffice: The second moment 
En>o ^"lt 11 ) < M for some constant M and, since the 
largest degree is 0(L 1 ^ b ~ 1 " > ) with high probability, 



C 

n(n — 2)—? — n(n — 2)v\{n) 



> e < (39) 



exp 



n 2J \MJ 



for n = 0(L 1 ^ b *)) as L — > oo. The magnitude of the 
error term in calculating 



^ n(n - 2)£>„/L 



(40) 



n>l 



turns out to be crucial. Writing the estimate (13) in 
a multiplicative form (this is, in fact, exactly the same 
bound as given in the appendix A of Ref. [26]) and using 
again the inequality e/(e) > e 2 /(2 + e) yields 



T 



> eu{(n) I < exp I —v\(n)L 



2 + e 



(41) 

so that we get a little sharper estimate than that given 
by the formula (14): 



C n = u[(n)L + ot^L 1 ' 2 ^). (42) 

Since v\(n) <v\(n)jr <~ rT h for large n, the error in the 
quantity (40) is 



0(L 



n 2 ^{n)L- 1 / 2 + 5 = o( K Li^h +s y (43) 



n=l 



Thus > 3.51 is enough to get rid of it. 



A. Phase diagram and critical exponents 

The simple approximation by independent degrees fails 
to give positive results even for the existence of the giant 
component for all the values of the parameter b. Now, 
as we start exploring the phase diagram, we present the 
uncertain cases as conjectures — they will be given some 
support by numerics exactly at criticality and analysis 
just below it, suggesting a cross-over between the sub- 
critical and non-trivial behaviors. The curves A-F of 
Fig. 2 show the critical curves for a set of values of 6, 
each representing a type of diagram different from the 
others. They will be commented one at a time as we 
discuss the corresponding values b. 

A remark on the calculations: Since many of the results 
discuss the behavior of the system in the high density 
limit, in that (f> < 1, the linear transformation formula 




FIG. 2: The curves A-F show phase diagrams for the Evans 
interaction with various values of b: From A to F, b — 
0, 3/2, 7/3, 8/3, 7/2, 9/2 and 8. The dotted line P is the curve 
S(p) of limiting points for 3 < b < 7 and the dotted line Q 
marks the value of density p, below which the birth of the 
giant component is possible only through zero-range conden- 
sation. 



for hypergeometric functions (15.3.6 of Ref. [25]) 

„, ,\ T(c)r(c - on - a 2 ) 

F(a 1 ,a 2 ;c;4>) = — r x 

F(c - ai)r(c- o 2 ) 

F(ai, a 2 ; ai + 02 - c + 1; 1 - (j>) 

r(c)r(ai + a 2 - c) 



(44) 



+ (1-0) 



c—ai —ci'2 ' 



r(oi)r(o 2 ) 

F(c — ai, c — a 2 ; c — a\ — a 2 + 1; 1 



0), 



valid for c / «i + Q2 + m, m e Z, with expansions of the 
hypergeometric functions on the right-hand side around 
0=1, proves to be very useful. It turns out that the 
transformation formula is not valid for integer b in our 
model, so that these situations must be dealt with sep- 
arately by repeatedly using the Gauss' relations for con- 
tiguous functions given in Ref. [25]. 

The 1 / p-behavior of r c corresponding to b = is plot- 
ted as the curve A of Fig. 2. This benchmark case has 
already been covered in section IV A, so we start explor- 
ing the phase diagram from non-zero values of b. There 
are no condensates for < b < 1, so that the results of 
the subcritical cases apply. Using the above transforma- 
tion formula and expanding in Eqs. (34) and (36) one 
gets for high densities that 



r c (p) 



1 



(2-b)p' 



(45) 



exactly the same as the prediction (32)! According to 
our calculations in section IV B this should be the only 
region, where the critical curve decays as a hrst inverse 
power of p. Indeed, at b = 1, 



z(0 = -(io g (i-0))M 



(46) 
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leading to a logarithmic correction to the first inverse 
power law. Furthermore, for 1 < b < 2, condensation 
does not occur, but we see a drastic change in the behav- 
ior of the critical curve, 



r c (p) 



1 



(6-i)r(2-6)r(6) V 

b\ p ) 



(47) 



in that it decays with a 6-dependent exponent greater 
than 1. The curve B in Fig. 2 with 6 = 3/2 is an ex- 
ample from this region. The case b = 2 involves again 
logarithmic corrections. We remark that high density 
approximations for the amplitude A of the giant com- 
ponent could have been computed from Eq. (19) in the 
same manner as for the critical curves. What is impor- 
tant is that, by the subcriticality in the zero-range sense, 
the growth of the giant component is linear in the vicin- 
ity of the phase transition, i.e. the critical exponent for 
the size of the giant component equals unity. 

There is a finite zero-range critical density p c = 1/(6 — 
2) for 2 < b < 3, but the second derivative of the parti- 
tion function Z(<f>), i.e. the second factorial moment of a 
^-distributed variable, still diverges in the limit (j> — > 1. 
This means that r c hits zero at p = 1/(6 — 2). Again from 
the expansions of the hypergeometric functions, we have 



r c (p) 



(b-2Y 



[(b-i)T(3-b)T(b)}— \b-2 



- P 



(48) 



The curves Ci (6=7/3) and C 2 (6=8/3) in Fig. 2 with 
zero-range critical points 3 and 3/2, and high density 
exponents 2 and 1/2 for the critical curves, respectively, 
are examples from this region. At the exact zero-range 
criticality, there is a giant component for any non-zero r, 
but the lack of proof prevents us from exactly calculating 
its size. Certainly something interesting is happening: 
Keeping 8 = r — r c (<j)) fixed and letting <f> —> 1, we have 



z(<i>)z'»(<i>y 



o, 

8(6-4) 
. 9(6-2)(6-3) 



for 6 < 4, 
5, for 6 > 4, 



(49) 



so that the linear growth is expected to be taken over 
by some other behavior. To set up a conjecture, we now 
calculate the leading order approximation as if the con- 
ditions of Molloy and Reed were satisfied. So we wish to 
analyze 



(3 = 1 



1 



f(2,2;2 + 6;(l-/3r)0) 
F(2,2;2 + 6;0) ' 

F(l,l;l + fe;(l-/?r)0) 
F(l,l;l + M) 



(50) 



(51) 



at <f> = 1. Expansions of the hypergeometric functions to 
the leading order in f3 then yield for the size of the giant 
component the expression 



1 



6-2 



r(6-2) 



r(2 - 6)r(6) 2 



(52) 



A o.oi 



o.ooi 




o.ooi 



r-r,. 



FIG. 3: Simulation results for the size of the giant component 
at the zero-range critical point for b = 2.25. The dash-dotted 
line is a numerical solution for the Eqs. (16) and (17). The 
dotted line shows the conjecture (52) for the leading order. 



as r — > 0, in that the critical exponent would be 1/(3 — 6). 
Results of Monte Carlo simulations, presented in Fig. 3 
for 6 = 2.25, indicate that the Eqs. (50) and (51) re- 
ally describe the size of the giant component correctly, 
although the convergence seems to be quite slow [27]. 

To give more support to our conjecture, we remark 
that it can be shown that the non-trivial term is already 
present in the expansion of A for <j> just below the ra- 
dius of convergence. This is obtained by expanding the 
hypergeometric function in Eq. (50) in such a way that 
1 — (3r<j) "C 1, in practice meaning that we are suffi- 
ciently far away from the critical curve measured in terms 
of 1 — <f), and the result is the same as Eq. (52) but with 
r multiplied by </>. Since the coefficient of the linear term 
vanishes in the limit, we then know that, at least for 
6 < 2.5, a crossover must exist at some small value of 
r that vanishes in the limit <p — > 1. The cases 6 > 2.5 
would require a higher order analysis of the size of the 
giant component below the condensation point in order 
to make sure that the coefficients of the analytical terms 
with powers between 1 and 1/(6—3) tend to zero as well. 

The critical curve for b = 3 is of the same type as for 
the 2 < 6 < 3 in the sense that it approaches zero as 
p tends to p c . The linear transformation formula being 
not valid for integer 6, one has to deal with the exact 
expressions for the derivatives of the partition function. 
This implies the conjecture 



A exp 

for the critical growth as r 
transition. 



1 

Ar 



(53) 



0, i.e. an infinite order 
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Now we enter the region 3 < b < 4, nearly half of 
which allows a rigorous analysis. Strictly speaking, we 
don't even know the exact position of the phase transi- 
tion for b < 3.5 at 4> = 1. This is because the limit of 
the critical curve r c as <p — > 1 is greater than zero by the 
existence of the second moment of a ^i-distributcd vari- 
able. More precisely, we have lim^i r c (</>) = (6 — 3)/4, 
this holding for all b > 3, so that, by p c = 1/(6 — 2), the 
end points define a curve 



S(p) = max < 0, 



1 (I 



4 \p 



1 



(54) 



on the (p, r)-plane. Below S(p), the transition curves 
are vertical and the giant component can be born only 
through condensation of edges. The way the curves r c (p) 
approach the zero-range critical point can again be eval- 
uated, 



6-3 



+ 



r(6) 2 r(4-6) 

16(6-3)T(6-3) 



(fr-2) 2 (b-3) 
36-5 



6-2 



(55) 

6-3 



The curve D of Fig. 2 shows an example with 6 = 7/2. In 
the figure, the curve P drawn with a dotted line is S(p), 
the curve of limiting points of r c . For 3.51 < b < 4 we 
can now prove that, at p = p c = 1/(6 — 2) as r — > r c (p) = 
(6 — 3)/4 from above, 



1 



16r(6- 2) 



6-2 V(&-3) 2 r(6) 2 r(2-6) 



6-3 



(56) 

thus giving the exponent 1/(6 — 3), exactly the inverse of 
the exponent for the critical curve. Remember that for 
2 < 6 < 3 we conjectured 1/(3-6) for A but (3-6)/(6-2) 
for r c (p), so that no such relation is valid for small values 
of 6. The case b = 4 again involves logarithms. 

We have two types of critical behavior to deal with left. 
The first one, occurring in the region 4 < b < 7, can be 
seen as a continuation to the preceding interval, since the 
cndpoints of r c lie on the curve S(p) with values strictly 
less than one because p c > 1/5. The change is in the 
exponents: Near the zero-range criticality we have 



rc{p) 



6-3 /9(6-3) \ (6-2) 2 (6-3) / 1 



- + 



4(6 - 4) 



-1 



36-5 



6-2 



P 



(57) 

in that r c forms a cusp with the vertical line at p c (see the 
curve E with 6 = 9/2 in Fig. 2). Moreover, we recover the 
linear growth of the giant component, that we predicted 
already from the subcritical analysis in formula (49), 



9(6-2)(6-3) 1 



(58) 



The region in the phase diagram that we have not yet 
discussed is b > 7. But this is trivial, since S(p) > 1 
for all values p < p c , and the phase transition curves 



are vertical lines positioned at the critical density of the 
zero-range process (curve F in Fig. 2). In other words, no 
giant component can exist without a condensate of edge 
ends on one of the vertices. 

We remark that the critical exponents for the growth 
of the giant component at the condensation point are ex- 
actly the same as obtained recently by Cohen et al. [28] 
for diluted scale-free networks. This is not a coincidence: 
In their model, from a random graph with a given power- 
law degree distribution, a fraction of the vertices with all 
the edges joined to them are removed. For the vertices 
retained in the graph this means that their degree dis- 
tribution is transformed exactly in the same manner as 
the distribution for the vertices in W in the operation 
of removing the dangling ends from our bipartite graph. 
The distribution v\ corresponding to the critical density 
being heavy-tailed, the similarity between the two prob- 
lems is obvious. In the dynamical model, however, we 
have to deal with the fluctuations and, as the analysis 
shows, they are not always easily controlled. 



VI. CONCLUSIONS 

We studied the existence and size of the giant com- 
ponent on one of the vertex sets in a dynamical model 
for bipartite multigraphs using the results of Molloy and 
Reed [12, 13]. Motivated by the concept of preferential 
attachment, we defined the movement of the edge ends on 
one side of the bipartite graph to be a zero-range process 
with jump rates bounded away from zero, in the station- 
ary state of which only weak correlations exist. A simple 
approximation by independent random variables allowed 
us to bound the fluctuations, and show that the degree 
sequences obtained this way are well-behaved in the sense 
of Molloy and Reed for zero-range processes with density 
below the critical value for condensation. As a bench- 
mark, it was shown that the case of non-interacting edge 
ends produced the well-known results of the G^M-model 
of random graphs. Generally, we find four types of crit- 
ical curves, depending on the asymptotics of the deriva- 
tives of the grand canonical partition function near the 
radius of convergence. A class of rate functions, for which 
the curves are of the same type as for the independent 
walkers, was also identified. The last sections of the ar- 
ticle were devoted to the Evans interaction, showing all 
the four types of phase diagrams as the parameter of the 
interaction is varied. Our simple use of the condition- 
ing device was strong enough to access some nontrivial 
parts of the parameter space but did not allow us to 
show that the conditions of Molloy and Reed are satis- 
fied for all values of the interaction parameter exactly 
at the zero-range condensation point. The critical expo- 
nents for the growth of the giant component were given 
as conjectures supported by Monte Carlo simulations in 
those cases. In the end, we discussed the connection with 
the model of diluted scale- free networks [28] . This implies 
that our model with the Evans interaction and exactly 
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at the zero-range condensation point belongs to the uni- 
versality class of percolation on graphs with power-law 
distributed degrees. The same set of exponents appears 
also in Ref. [29] , where a general, phenomenological Lan- 
dau theory for phase transitions on scale-free networks 
is constructed. The simulation results suggest that the 
simple Landau theory is not broken by the fluctuations 
in our model. 

Finally, we would like to mention a few points that 
would be interesting to analyze further. First of all, our 
simple approximations by independent variables failed to 
give estimates sharp enough to completely cover the parts 
of the parameter space with condensates for the Evans 
interaction. We believe that the problem is purely tech- 
nical and could be removed by more sophisticated use of 
the conditioning device. Secondly, no results 'inside' the 
phase transition, i.e. no further than o(l) from the criti- 
cal point, exist for graphs with a general degree sequence. 
Therefore, the behavior near r = r c remains to be dis- 
covered. What comes to the underlying particle system, 
there are still open questions even concerning the stat- 
ics of zero-range processes, especially for rate functions 
exhibiting slow decay [16]. Furthermore, only station- 
ary properties of the model were discussed in this article. 
Recently, advances have been made towards understand- 
ing the formation of condensates in zero-range processes 
as a function of time [15, 30], so that one is tempted to 
consider the temporal scalings for graph- valued processes 
as well. A grand canonical generalization of the present 
model that would allow for a variable number of external 
agents would also be of interest. 
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APPENDIX A: DEFINITIONS FROM [12, 13] 

We give the definitions of Molloy and Reed in our no- 
tation. An asymptotic degree sequence is a sequence of 
integer- valued functions V = Dq{L),D\{L), . . . such that 
D n {L) = 0, for n > L and £„> D n {L) = L. 

Given an asymptotic degree sequence V, T>l is set to 
be the degree sequence ci, C2, . . . , Cl, where Cj > cj +1 
and \j : Cj = n\ = D n (L). Let CI-d l be the set of graphs 
with vertex set 1, . . . , L with degree sequence T>l. 

An asymptotic degree sequence V is well-behaved if: 

1. T> is feasible and smooth, i.e. Qx> L ^ and there 
are constants A„ such that Wuil^oo D n (L)/L = X n . 

2. n(n — 2)D n {L) / L tends uniformly to n(n — 2)A„, 
i.e. for all e > there exists L e such that for all 



L > L e and for all n > 

,n(n - 2)D n (L) , , 
P ^ " V > - n(n - 2)A„| < e 

3. 

C(V)= lim V n(n-2)D n (L)/L 

n>l 

exists , and the sum approaches the limit uniformly, 
i.e.: 

(a) If C(V) is finite then for all e > 0, there exist 
n* , L e such that for all L > L e : 

^(n-2)D n (L) _ cm<e 

71=1 

(b) If C(D) is infinite then for all T > , there 
exist n*, L e such that for all L > L e : 

| ^ n(n - 2)D n {L) | ^ y 

n=l ^ 



Furthermore, V is called sparse if nD n (L)/L = K + 
o(l) for some constant K. 

In the theorems of Molloy and Reed, T> is assumed 
to be a well-behaved sparse asymptotic degree sequence 
with the property that there is an e > such that for 
all L and n > L 1 / 4-6 , D n = 0. In proving some state- 
ments concerning the cycles 1/4 — e has to be replaced 
by more restrictive 1/8 — e. The figure 1/4 is needed for 
the multigraph of a random configuration to be simple, 
and is therefore not required in our analysis. In fact only 
1/2 — e, which is high enough for our purposes, is needed 
for the statements about the existence and size of the 
giant component to hold. 

Notice that our sequence D defined in section III can- 
not satisfy the conditions almost surely (i.e. with prob- 
ability equal to 1) but with high probability only. This 
does not matter, since all the results are given in the 
probabilistic setting anyway, and basically means a one 
more choice of L e in the proofs. 
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